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PACS 05 . 45 . Xt - Synchronization; coupled oscillators 
PACS 89 . 75 . He - Networks and genealogical trees 

Abstract - We introduce a model to study the effect of degree-frequency correlations on syn- 
chronization in networks of coupled oscillators. Analyzing this model, we find several remarkable 
characteristics. We find a stationary synchronized state that is (i) universal, i.e., the degree of 
synchrony, as measured by a global order parameter, is independent of network topology, and (ii) 
fully phase-locked, i.e., all oscillators become simultaneously phase-locked despite having different 
natural frequencies. This state separates qualitatively different behaviors for two other classes of 
correlations where, respectively, slow and fast oscillators can remain unsynchronized. We close by 
presenting analysis of the dynamics under arbitrary degree-frequency correlations. 
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Introduction. — The research of emergent collective 
behavior in large ensembles of interacting dynamical sys- 
tems represents a large and important area of complexity 
theory [lj^j. Studying synchronization of coupled oscil- 
lators has proven to be particularly useful in modeling 
complex systems and uncovering generic mechanisms be- 
hind synchronization processes. Examples include simul- 
taneous flashing of fireflies [5] , cardiac pacemaker cells [6] , 
circadian rhythms of mammals |7j, collective oscillations 
of pedestrian bridges j8j, and chemical oscillators [9 . In 
many cases, the interactions between oscillators can be 
described by a complex network. To gain insight into the 
mechanism behind synchronization, Kuramoto proposed 
to model the state of each oscillator n by a phase variable 
n 10 . When placed on a network, the dynamics of n is 



governed by 



N 



K ^2 A nm Sin(0 m -d n ), 



(1) 



where u> n represents the natural frequency of oscillator 
n, K is the global coupling strength, and [A nm ] is the 
adjacency matrix that encodes the network topology of 
the underlying system (n, m = 1, 2, . . . , N). 

Although network topology plays a vital role in deter- 
mining synchronization 11-20 , the question of how it in- 



fluences synchrony is not completely understood. In re- 
cent years researchers have started to explore the effect 
of correlations between oscillator frequency uj n and de- 
gree k n — Yl m =l ^-nm and observed that in some cases, 
such correlations can give rise to enhanced synchroniz- 
ability 11 and the emergence of explosive synchronization 
What, then, is the effect of degree-frequency 
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events 

correlations on synchronization in general? 

We address this question by analytically and numer- 
ically studying synchronization in undirected networks 
(i.e., those for which A nm = A mn ) with general degree- 
frequency correlations. These correlations may be char- 
acterized by the joint probability distribution of degrees 
and frequencies P(k,w), which we assume to be symmet- 
ric about u) = 0, i.e., P(k,— u) — P(k,w). In the classi- 
cal (uncorrelated) network Kuramoto model 



13 15 , the 



frequencies and degrees are chosen independently, so that 
the distribution P{k, w) can be written as a product of the 
frequency and degree distributions, P(k,oj) = P(k)g(w). 
In this Letter we propose a framework to study synchro- 
nization in the general case and present detailed results 
for the case in which the joint distribution is given by 
P(k,u) = P(k)[6(u - ak fj ) + 6(lu + ak f) )]/2, i.e., 



±akt 



(2) 



where a, (3 characterize the correlation and the positive 
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and negative signs are chosen with equal probability to 
maintain zero mean frequency as N — > oo. This particular 
form of P(k, u>) is chosen as an illustrative example and is 
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closely related to a model studied numerically in Ref. 
(see their footnote [24]). 

This simple model can be used for analyzing the influ- 
ence of degree-frequency correlations on the synchroniza- 
tion of coupled oscillators and exhibits rich dynamics. We 
note that a can be scaled out of eqs. ([T]) and ^ by letting 
t h-> t/a and K i-> aK. We will therefore use a = 1 in all 
figures presented in this Letter. Thus, the free parameters 
are the coupling strength K, the adjacency matrix L4 nm ], 
and the correlation exponent f3. We will consider positive 
correlations and refer to (3 — 1, (3 < 1, and j3 > 1 as linear, 
sub-linear, and super-linear correlations, respectively. 

To measure the degree of synchrony, we introduce the 
following order parameters. The local order parameter r n 
for oscillator n, which quantifies the degree of synchrony 
among the neighbors of node n, is defined by r n e l ^ n = 
J2 m Anme iem , where ip n is the local mean phase. The 
global order parameter is defined by R = N^ 1 J2 n JT anc ^ 
measures the degree of synchrony over the entire network. 

Description of solutions. — We now briefly describe 
the dynamics of the steady-state behavior of the system 
defined by eqs. and We begin by describing how 
synchronization changes as the coupling strength K is var- 
ied. In fig. [I] we plot data from simulations on an Erdos- 
Renyi (ER) network 21 of size N = 1000 with link prob- 
ability p = 0.1, using a correlation exponent (3 = 1. Fig- 
ure[l|a) shows that as the coupling strength K increases, 
the time-averaged order parameter R also increases to- 
wards the value of 1, as expected. Notably, this R-K 
curve exhibits two transitions, one at the critical coupling 
strength K — K\ ss 0.2, and the other one at the critical 
coupling strength K = K 2 ~ 2 (indicated with vertical 
dotted lines). This is in sharp contrast to the usual R-K 
curves where a single transition is observed [l3] . 

As shown in fig. [lja), these two critical coupling 
strengths separate three regimes which we denote as inco- 
herent (I), standing wave (SW), and stationary synchro- 
nized (SS) states. For K < K\, R ps and the system 
is incoherent, consisting of oscillators that drift indepen- 
dently. For Ki < K < K 2 , the network exhibits S W solu- 
tions characterized by the emergence of two synchronized 
clusters traveling with opposite angular velocities. Such 
SW solutions result in the oscillating behavior of R(t), 
shown in fig. [ljb). The distribution of phases p(9) corre- 
sponding to the maximal and minimal R(t) values [e.g., 
as indicated by the green circle and red cross in fig. [ljb) , 
respectively, and shown as dashed lines in fig. 2(a)] are 
depicted in fig. [ljc) . Note that R(t) achieves its maxi- 
mum when the distributions of phases for the two clusters 
overlap (dashed green) and achieves its minimum when 
they lie on opposite sides of the unit circle (dot-dashed 
red). Finally, for K > K 2 an SS state emerges, yielding a 
time-invariant R(t) ~ 1 [fig. Qd)]. 
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Fig. 1: (Colour online) Transition from incoherence to coher- 
ence under the presence of degree-frequency correlation. All for 
an ER network with parameters N = 1000, p = 0.1, and /3 = 
1. (a) Time-averaged (solid blue) and minimum/ maximum 
(dashed green) R versus K. (b) Time series R(t) for an SW so- 
lution using K = 1.9. (c) Distribution of phases p(6) of phases 
at maximum (dashed green) and minimum (dot-dashed red) 
R(t) values [times denoted by the green circle and red cross in 
(b)]. (d) Time series R(t) for an SS solution using K = 2.1. 



The SS state exhibits remarkable characteristics. In 
particular, as we will see, with a linear correlation the 
critical coupling strength for the onset of global synchro- 
nization is K~2 = 2a, a universal value that is independent 
of detailed network topology. The steady-state degree of 
global synchrony R as a function of K also turns out to 
be universal in the case of a linear correlation. In sharp 
contrast, when there is no degree-frequency correlation or 
when such a correlation is nonlinear, network structure 
plays a vital role in determining both the critical coupling 

Fur- 



strength and degree of global synchrony R 13 15 
thermore, in the absence of a degree-frequency correlation, 
only a fraction of the oscillators become phase-locked. The 
oscillators that are not phase-locked drift indefinitely and 
typically have either low degrees or high frequencies [13] . 
However, for a linear degree-frequency correlation (/3 = 1), 
whenever the system exhibits global synchrony, all oscil- 
lators are locked, which we refer to as full phase-locking. 
For nonlinear correlations, we find (through both analyti- 
cal and numerical approaches) that when the correlation is 
super-linear (sub- linear) , drifting oscillators typically exist 
and are those with high (low) degrees. A linear correlation 
thus represents a perfect balance between each oscillator's 
topological (degree) and dynamical (frequency) proper- 
ties. We illustrate this in fig. [2j where we show locked 
(blue) and drifting (yellow) oscillators from a real simula- 
tion of a network of size N = 16 for sub-linear, linear, and 
super-linear correlations (left to right). Note that for the 
sub-linear correlation only oscillators with small degrees 
(k n = 2) drift, while for the super-linear correlation only 
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Fig. 2: (Colour online) Illustration of phase-locking for sub- 
linear, linear, and super-linear correlations in a network of 
size N = 16. Circle radii are proportional to degrees with 
locked and drifting oscillators colored blue and yellow, respec- 
tively. Simulation parameter values are j3 = 0.8, 1, 1.2 and 
K = 1.35,2.1,2.5, respectively. 



oscillators with large degrees (k n > 8) drift. The case of 
a linear correlation corresponds to full phase-locking. 

Standing wave solution. — The existence of the SW 
state can be understood by noting that the frequency dis- 
tribution of the oscillators is bimodal, a property that has 
been previously shown to produce SW states for systems 
lacking degree-frequency correlations 22 . For example, 
when a, ft = 1 in eq. ([2]), the frequency distribution g(uj) 



is simply the mirror-reflected version of the degree dis- 
tribution, g(uj) = [P(— (J) + P(lo)]/2. Thus, a unimodal 
P(k) [in the case of an Erdos-Renyi network, peaked at 
k = p(N— 1)] naturally gives rise to a bimodal g(uj), which 
is expected to lead to an SW solution when the separation 
between the two peaks of g(oj) is large enough compared 
to the width of the distribution [22] . 

To begin the analysis of the SW solution, we will an- 
alyze separately the degree of synchrony in the clusters 
of oscillators with positive and negative frequencies. To 
this end, we introduce positive/negative local and global 
order parameters r^e 1 ^™ = ^2 Um ^ A nm e l6m and BF = 
N~ x J2 n r n An; where k„ is the sum of link strengths con- 
necting oscillator n to oscillators with positive/negative 



frequencies, fc± = E^o 
order parameters, eq. (fl| < 

6 n = uj n + K[r+ sin(^>+ 



A nm . Using the modified local 
can be rewritten as 



sin(V>„ 



.)]■ (3) 



We now assume that synchronized oscillators are divided, 
according to the sign of their frequency uj n , into two clus- 
ters that rotate in opposite directions with angular veloc- 
ity ±f2, so that — ±Qt. Assuming u> n > and moving 
to a rotating frame of coordinates, we define cf> n = 9 n — fit, 
and obtain 

K = (w„ - O) - Kr+ sin(0„) - Kr~ sin(0„ + 2fii). (4) 

For not too small, the last term in this equation oscil- 
lates rapidly around zero compared to the first two terms 
and can therefore be approximately averaged out. (Later 



we will discuss when the value of f2 we find in our analysis 
is consistent with this assumption.) We will now look for a 
solution in which the values of the local order parameters 
r+ are approximately time-independent. We note that this 
occurs when oscillator degrees fc+ ss k n /2 are large enough 
that fluctuations may be neglected (see [13] for a discus- 
sion) . Accordingly, we neglect the last term in eq. Q , take 
r+ to be independent of time, and find that oscillator n 
locks with the positive cluster if \u) n — f2| < Kr£, in which 
case we have that sin(</>„) = ; otherwise, it drifts 

indefinitely. Due to the symmetry of the frequency distri- 
bution, drifting oscillators (as a whole) do not contribute 
to the degree of local or global synchrony [13] , allowing us 
to rewrite the local order parameter as 



E 

U! n >0, 



(5) 



Now, since exactly fc+ terms contribute to the order pa- 
rameter r+, we propose that r+ is proportional to fc+. 
This approximation has been validated numerically for 
this and other network-coupled oscillator systems, but is 
expected to break down for small r+ in very heterogeneous 
networks, e.g., networks with a scal e-free ( SF) degree dis- 
tribution P(k) oc fc -7 with 7 < 2.5 
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Therefore, we 

expect the following theory to be valid only for relatively 
homogeneous networks. Given the definition of R + , we 
set r+ = i? + fc+. Recalling that u> n = ak^ for uj n > 0, we 
separate eq. ^ into its real and imaginary part to obtain 
self-consistent expressions for R + and fi, 



I? 4 



(k)- 



N 



2\ak p „ 



-Q.\<KR+k m 



4(aC - ny 



{KR+k r , 



n = 



^2\ak?, 

! ^ 



-n\<KR+k r 
-Q\<KR+k„ 



k'L 



(6) 
(7) 



where (k) — ^ k n /N and we have also used fc+ w k m /2. 
For large N, eqs. ^ and ^ can be approximated by 



R+ = (k)- 1 



2\akf-a\<KR+k 



/ 4(afc' 9 - m 2 „ 



f2\akf>- 



n\<KR+ 



k P(k)k?dk 



^2\a 



\akP-n\<KR.+ 



k P(k)dk 



(8) 
(9) 



A similar argument would show that R~ satisfies eq. 
Eqs. ^ and ^ give the degree of synchrony in each clus- 
ter and must be solved self-consistently. In general eqs. ^ 
and ([9| need to be solved numerically. 

In the case of j3 — 1, it is possible to find analytically 
the critical value K\ corresponding to the onset of the SW 
solution. To do this, we substitute z = 2(ak — il)/KR + k 
in eq. (18]) and let R + — > + , obtaining a critical coupling 
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Fig. 3: (Colour online) Degree of synchrony within positive 
clusters R + versus coupling strength K for an ER network 
with p = 0.1 (blue circles) and a SF network with 7 = 5.0 
and ko = 50 (red crosses), both of size N = 1000. Theoretical 
predictions for R + given by eqs. ^ and ^ are plotted in 
dashed black. Critical coupling strengths Ki are marked with 
vertical dotted lines. 



strength of K\ — 4a 3 (fc) /7rQfP(Oi/a), where f2i is the 
group angular velocity at onset. If P(k) is unimodal and 
has a peak at an intermediate k value, e.g. for an ER 
network, then expanding eq. |9| about R + = yields 
the condition P'(Q,\/a) — 0. For an ER network with 
mean degree (k) = (N — l)p, this yields f2i = a(k), 
K\ — 4a/ir(k)P((k)). For a monotonically-decreasing 
distribution P{k) with minimum degree k , e.g., a SF 
network with minimum degree ko, it can be shown that 
Oi = afco, which yields a critical coupling strength of 
K\ = 4a(k) /TrkyP(k()). We note that, at onset, the period 
of oscillation of the last term in eq. Q is n/Qi. On the 
other hand, the timescale of evolution associated with the 
first two terms is 27r/(oj n — f2i). Therefore, to neglect the 
last term in eq. Q we require 2Qi 3> u) n — fl±. For a dis- 
tribution peaked at k = k we require, using f2j = ak and 
Lj n = ak n , that 2k k n — k. Therefore we strictly require 
2k 3> max„(fc„ — k). A somewhat less restrictive require- 
ment, which guarantees the condition is valid for most of 
the oscillators, is 2k 3> rms(fc„ — k). In any case, our the- 
ory for the onset of the standing wave solution is restricted 
to networks with a homogeneous degree distribution (e.g. 
not SF networks). 

We numerically verify these results by simulating 
eqs. ([I]) and Q with (3 — 1 over a range of K for an 
ER network with p = 0.1 and a SF network with 7 = 5 
and fco = 50 (all SF networks we use in this Letter were 
generated using the configuration model 23 ). Both net- 
works are of size N = 1000. Resulting R + for the ER and 
SF networks are plotted in blue circles and red crosses, re- 
spectively, in fig. [3] We note that corresponding R~ values 
were indistinguishable from R + . Theoretical predictions 
obtained by solving eqs. Q and ^ are plotted as dashed 
black curves. Critical values Ki for each network are indi- 
cated by vertical dotted lines. Results from simulations on 
the ER network are predicted well by our theory. While 
our theory is not expected to apply to SF networks (see 
previous paragraph) we find reasonable agreement for the 
SF network with 7 = 5. The agreement does break down 



Stationary synchronized solution. We now 

present an analysis of the SS solution. Using the definition 
of the local order parameters we rewrite eq. ([!]) as 



u> n + Kr n sm(ip n 



.)• 



(10) 



We now look for solutions where (i) the synchronized clus- 
ter has zero mean frequency and (ii) local order param- 
eters r n are approximately time- invariant. Oscillator n 
then becomes phase-locked if \ui n \ < Kr n , in which case 
sin(0 n — ip n ) = u> n /Kr n ; otherwise it drifts indefinitely. 
Due to the symmetry of the frequency distribution, drift- 
ing oscillators (as a whole) do not contribute to the degree 
of local or global synchrony [l3j , allowing us to rewrite the 
local order parameter as 



%\<Kr„ 



(11) 



We now look for solutions that satisfy the following con- 
ditions. First, following 13 we set tp n = ip m for all n,m. 



Second, as in the analysis of the SW solution, since exactly 
k n terms contribute to the order parameter r n , we propose 
that r„ is proportional to the degree k n , i.e. r„ = Rk n . 
We note that this holds extremely well even for very het- 
erogeneous networks because for SS solutions r n /k n 1. 



Under these two assumptions, eq. (11) becomes 



Rk n — ' A nm \ 

\u m \<KRk m 



1- 



\KRk n 



(12) 



For the linear correlation ((3 = 1) the dependence of 
both the summation condition and the square-root term 
on k m (and u> m ) vanishes. Looking for the synchronized 
state, we sum eq. (12) over all nodes and find that, after 
some simplification, 



R 



\ 



1± ^ 1 



i K 2 



(13) 



where the + (— ) sign represents a stable (unstable) so- 
lution (numerically determined). This branch of station- 
ary synchronized solutions appears at Ki = 2a in the 



form of a saddle-node bifurcation. Note that in eq. (12), 



since the square-root term becomes constant, the remain- 
ing J2 m Anm term, which encodes the network topology, 
reduces to the degree k n which is balanced by the left side 
of eq. (12). Thus, the degree of global synchrony given 
by eq. ( 13 ) and the critical coupling constant K2 — 2a 



at which the SS solution appears are independent of the 
detailed structure of the network, which we refer to as uni- 
versality. This surprising result is found to hold even for 
networks with degree-degree correlations. 

We numerically verify these results by simulating 
eqs. |T]) and Q with ft = 1 over a range of K for an ER 
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Fig. 4: (Colour online) Degree of global synchrony R versus 
coupling strength K. (a) Several networks with linear corre- 
lations, P = 1. Networks used are ER with p — 0.1 (blue 
circles), and SF with 7 = 2.5 (red crosses) and 3.5 (green tri- 
angles), both with fco = 10, all of size TV = 1000. Theoretical 
prediction given by eq. ( 13 1 in dashed black, (b) Nonlinear 
correlations /3 = 0.9 (blue circles) and /3 = 1.1 (red crosses) 
on a SF network with 7 = 3, fco = 50. Theoretical predictions 



given by eq. (151 in dashed black. 



network with p = 0.1 and two SF networks with 7 = 2.5 
and 3.5 and fco = 10. All networks are of size N = 1000. 
Resulting R values for the ER network and SF networks 
with 7 = 2.5 and 3.5 are plotted in blue circles, red crosses, 
and green triangles, respectively, in fig. [4^a) . The theoreti- 
cal prediction given by eq. ( 13 1 is plotted in dashed black. 
The critical coupling strength K 2 — 2a is indicated by 
the vertical dotted line. Results from simulations are pre- 
dicted very well by our theory, confirming that the detailed 
network topology is not necessary to describe K 2 and R 
in the SS state for linear correlations. We note that, as 
opposed to our theory for the SW solution, here we do not 
need to assume a homogeneous degree distribution. 

For nonlinear correlations {ft 7^ 1), eqs. ^ and (12) 
yield, after summing over n, 



R = 



(k)' 



N 



E 



ak%.<KRk„ 



where (fc) = J2"=i k n/N. For large N, we can approxi- 
mate eq. ([l4| with the integral 



1 - 



f-km 



KRk„ 



(14) 



R=(ky 



akP<KRk 



P(k)k\ 1 - 



akf 3 
KRk 



dk. (15) 



In general, eq. ( 15 ) needs to be solved numerically. 



The critical coupling strength K 2 where the station- 
ary synchronized solution is born can be found by solving 



eq. ( 15 ) for the minimum K value where R > is a solu- 



universal value independent of the network topology. In- 
creasing (decreasing) j3 effectively spreads (contracts) the 
set of natural frequencies, therefore impeding (promoting) 
synchrony and increasing (decreasing) K 2 . 

We numerically verify these results by simulating 
eqs. ([TJ and ^ with /3 7^ 1 on a SF network with 7 = 3 
and minimum degree fco = 50. Resulting R for j3 = 0.9 
and 1.1 are plotted as blue circles and red crosses, re- 
spectively, in fig. |4^b). Theoretical predictions for R and 
the critical coupling strength K 2 , both obtained by solv- 



ing eq. (151, are plotted as dashed black and vertical dot- 



ted black curves. Results from simulations are predicted 
very well by our theory. For networks which violate our 
assumptions by having smaller minimum degrees, e.g., 
fco = 10, we found that K 2 as observed from simulations 
is slightly smaller (larger) for f3 < 1 ((3 > 1) than those 



predicted by eq. (15) (simulations not shown). 



To further explore the dependence of K 2 on network 
characteristics, we consider SF networks and numerically 



solve eq. ( 15 ) to find K 2 given a correlation exponent j3 



and degree exponent 7. Setting the minimum degree fco = 
50, we plot K 2 as a function of 7 in fig. [5] for increasing 
values of /3 € [0.8, 1.2], from bottom to top. We see that 
for j3 < 1, we have K 2 < 2a, and for j3 > 1, we have 
K 2 > 2a. As the networks become more heterogeneous 
(i.e., 7 decreases) K 2 curves upward (downward) for f3 > 1 
(/? < 1), while K 2 = 2a remains constant for /3 = 1. 

Having analyzed the SS state, we finally revisit the novel 
phase-locking behavior introduced in fig. [2] Recall our ob- 
servation that the linear (j3 = 1) correlation produces full 
phase-locking, implying that there are no drifting oscilla- 
tors. In fact this was observed to be a critical case separat- 
ing the contrasting phase-locking behaviors of sub-linear 
and super-linear correlations, for which there exist drifting 
oscillators with low and high degrees, respectively. This 
interesting phenomenon can be explained by the locking 
criterion afc^" 1 < KR in eq. (14), assuming that K > K 2 . 



For super-linear correlations (/3 > 1), oscillators with de- 
gree fc < {^£')^ T become locked, while oscillators with 
high degree and frequency drift, a scenario similar to what 
has been observed in previous work 13 . For sub-linear 



correlations (/3 < 1), the phase- locked population consists 
of oscillators with degree fc > (t^) 1 ^, thus leaving os- 
cillators with low degree and frequency drifting. These 
two qualitatively different behaviors are separated by the 
critical case of linear correlations (/? = 1) for which the 
dependence on fc disappears and the oscillators either all 
drift or all phase-lock. 

General Correlations. — We finalize our analysis 
by noting that, although in this Letter we focused on a 
specific form for the degree- frequency correlations [i.e., 
eq. Q], in the general case of a joint distribution P(fc, oS) 
symmetric about U) = 0, our analysis still holds and re- 



tion. Recall that for f3 — 1 we have K 2 — 2a, which is a suits generalize. For the SS solution, we find that eq. (15) 
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with degree-frequency correlations. This framework may 
potentially be used to optimize the synchronization prop- 
erties of networks, which have been recently realized ex- 
perimentally [25| . 
A recent preprint 



26 independently studied additional 



aspects of correlations between degrees and frequencies. 



We note that eq. ( 16 1 is complementary to these results. 



Fig. 5: (Colour online) Critical coupling strength K2 obtained 
from eq. ( |15[ ) as a function of 7 for SF networks with k — 50 
for several values of j3. 



is replaced with 



R=(ky 1 f f P(k,w)k x 1 

JO J\u>\<KRk y 



(KRk) 



rdiudk. 



(16) 

For a general distribution P(k,uj), the SW solution will 
not appear if the distribution of frequencies is not suffi- 
ciently bimodal. If the distribution of frequencies is suffi- 
ciently bimodal, then we may replace eq. (JsJ) with 



I? 4 



2\ui-Q\<KR+k 



P(k,Lo)k 



'l - ^ -4wdu:dk. 

(KR+k) 2 



(17) 



Conclusion. — In many applications of network- 
coupled dynamical systems, a central questions is how 
the dynamics and network structure interact and give rise 
to emergent collective behavior 11-20 . For instance, in 



many systems the contribution of the network structure 
is encapsulated in one or more eigenvalues and eigenvec- 
tors of the network adjacency [l6j or Laplacian matri- 
Here we find that if the degree-frequency cor- 
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ces 

relation is chosen appropriately, then the network struc- 
ture has virtually no influence on the synchronization 
properties of the resulting dynamics. 

Full phase- locking, i.e., the simultaneous entrainment of 
all oscillators, in heterogeneous oscillator systems is also 
a novel finding. Typically, an extremely large value of 
K is needed to entrain all the oscillators in a large net- 
work when the oscillators are heterogeneous 13 , 24 . How- 
ever, in the presence of a linear degree-frequency correla- 
tion, all oscillators become phase-locked simultaneously as 
the coupling constant passes the critical value for global 
synchrony, K-z- This unexpected phenomenon emerges 
despite the presence of strong heterogeneity in both the 
network structure and oscillator dynamics. Another re- 
markable observation is that, for sublinear correlations, 
the locked oscillators are those with a frequency which is 
most different from the mean. 

In addition to the analysis of the particular case in 
eq. ([2]), we have presented a general formalism to an- 
alyze the synchronization of network-coupled oscillators 
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